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Abstract 

We report theoretical and numerical evaluations of the phase diagram for a model of patchy 
particles. Specifically we study hard-spheres whose surface is decorated by a small number / of 
identical sites ("sticky spots") interacting via a short-range square- well attraction. We theoreti- 
cally evaluate, solving the Wertheim theory, the location of the critical point and the gas-liquid 
coexistence line for several values of / and compare them to results of Gibbs and Grand Canonical 
Monte Carlo simulations. We study both ordered and disordered arrangements of the sites on the 
hard-sphere surface and confirm that patchiness has a strong effect on the phase diagram: the gas- 
liquid coexistence region in the temperature-density plane is significantly reduced as / decreases. 
We also theoretically evaluate the locus of specific heat maxima and the percolation line. 
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I. INTRODUCTION 



Patchy particles are particles interacting via a limited number of directional interactions. 
The anisotropy of the interaction leads to collective behaviors different from those of simple 
liquids. Gelation^^^, gas-liquid phase separation^, crystallization^^ and clustering are 
strongly affected by patchiness^^kLiLL^yi. Recently, a new generation of colloidal particles 
with chemically or physically patterned surfaces has been designed and synthesized in the 
attempt to provide valence to colloids^^^yyy&^Si. This relevant synthesis effort aims to 
generate super atoms — atoms at the nano and micro-scopic level — in order to reproduce 
and extend the atomic and molecular behavior on larger length scale. It also offers the 
possibility to export the supra-molecular chemistry ideas^i^ 1 ^ to new colloidal materials, 
opening the new field of supra-particle colloidal physics. Thus, a general effort to develop 
a deeper understanding of self-assembly and to construct a more unified theoretical under- 
pinning for this technologically and scientifically important field is crucial. The outcome of 
this effort may also have an impact in our understanding of the phase behavior of protein 
solutions, due to their intrinsic patchy character— i2£i2&21. 

Our recent work- has shown that the Wertheim theory^^, describes rather well the crit- 
ical properties of particles decorated on their surface by a predefined number of attractive 
sites. The Wertheim theory is a thermodynamic perturbation theory introduced to describe 
association under the hypothesis of a single-bond per patch, which means that an attractive 
site on a particle cannot bind simultaneously to two (or more) sites on another particle. The 
single-bond per patch condition can be naturally implemented in colloids by choosing an 
appropriate small ratio between the range of the attractive patches and the particles size. 
The single-bonding condition results also from DNA complementarity^ 1 ^ as well as from 
complementary "lock-and-key" interactions associated to biological specificity^ 1 ^. These 
types of interactions provide a versatile way of controlling inter-particle binding. An exten- 
sion of the Wertheim thermodynamic perturbation theory to interpret and/or predict the 
behavior of a wide range of substances with potential industrial applications is provided by 
the statistical associating fluid theory (SAFT)^^ and by its developments^ 1 ^. 

In previous works we have shown^ 1 ^ that for patchy colloidal particles with a small 
number of sticky sites the critical point of the gas-liquid phase separation moves towards 
small packing fraction (</>) and temperature (T) with decreasing the number of patches. 
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According to this study, liquid phases of vanishing density can be generated once a small 
fraction of polyfunctional particles is added to a system of bifunctional ones. Indeed, the 
study of binary mixtures of patchy particles with different functionalities allows to explore 
also the range of non-integer valence down to 2. This means that with the new generation of 
non-spherical sticky colloids^ 1 ^, it should be possible to realize "empty liquids"-, i.e. states 
with an arbitrarily small occupied packing fraction at temperature lower than the liquid- 
gas critical temperature. The shift with valence of the critical point, both in density and 
temperature, leads to substantial changes in phase behavior with branching: the reduction 
of the number of bonded nearest neighbours is accompanied by an enlargement of the region 
of stability of the liquid phase in the (T, 0) plane. This fact could favor the establishment, 
at low T and at small <ft, of homogeneous disordered materials, i.e. equilibrium disordered 
states in which particles are interconnected in a persistent network. At such low T, the 
bond-lifetime will become comparable to the experimental observation time. Under these 
conditions, it should be possible to approach dynamical arrested states continuously from 
equilibrium and to generate a state of matter as close as possible to an ideal gel^. 

In this article we extend the preliminary study of Ref.- reporting a Monte Carlo inves- 
tigation of the /-dependence of the critical point location for a model with a disordered 
arrangement of patches. The present study confirms the trend discussed in Ref.- for the 
corresponding ordered case. To evaluate the role of the valence on the coexistence region, 
we also numerically investigate the shape of the gas-liquid binodal line for the ordered case 
and compare it with theoretical predictions based on the Wertheim theory. Finally, we 
analytically calculate several equilibrium properties, such as the energy per particle, the 
specific heat, the extent of polymerization and the percolation line, to get insights on their 
/-dependence. 

We find that the reduction of valence is accompanied by a significant shift of the coex- 
istence curve towards low temperature and density. The percolation line is always found 
to lie above the critical point, merging with the gas-liquid spinodal at low density p. The 
liquid state is thus always characterized by an infinite spanning network. This confirms the 
possibility of observing, for large attraction strengths, dynamical arrested states driven by 
bonding (as opposed to packing) in single phase conditions, i.e. homogenous arrested states 
at low density. 

We also provide in the Appendix A a physical insight of the Wertheim theory by showing 
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that the theoretical expression for the free energy in Ref.— <22 is formally equivalent to the free 
energy of a system of non-interacting clusters distributed according to the Flory-Stockmayer 
cluster size distribution 39 . 

II. THE MODEL 

We focus on a system of particles modeled as hard-spheres of diameter a, whose surface 
is decorated by / bonding sites at fixed locations. Sites on different particles interact via a 
square-well potential. The interaction V(l,2) between particles 1 and 2 is 

/ / 

V(l, 2) = V HS (r 12 ) + V sw(r%) (1) 

i=i j=i 

where the individual sites are denoted by i and j, Vhs is the hard-sphere potential, Vgw{ x ) 
is a square-well interaction (of depth — uq for x < 5, otherwise) and ri2 and are 
respectively the vectors joining the particle-particle and the site-site (on different particles) 
centers. Geometric considerations for a three touching spheres configuration show that the 
choice 5 = 0. 5(^5-2^-1 )er ~ 0.119cr guarantees that each site is engaged at most in 
one bond. Hence, with this choice of 5, each particle can form only up to / bonds. We note 
that in this model bonding is properly defined: two particles are bonded when their pair 
interaction energy is -uq. Distances are measured in units of a. Temperature is measured in 
units of the potential depth (i.e. Boltzmann constant k% = 1). 

We study two different arrangements of the / sites on the particles surface. In the first 
case sites are arranged in a regular structure (see Fig. 1 of Ref.-). In the second case, the 
distribution of the sites is random and different for each particle. In this latter case, the 
only constraint on the site position is formulated on the basis of a minimum distance d m i n 
criterion between different sites on the same particle: the choice of d m i n aims to minimize 
the possibility of double bonding between the same pair of particles as well as the shading 
of a bonding site by the presence of a nearby bonded site. We choose d min = 0.4. 

III. THE THEORY 

The first-order thermodynamic perturbation Wertheim theory^^ 1 ^ provides an expres- 
sion for the free energy of particles with a number / of attractive sticky sites on their 
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surface, independently from the specific geometric arrangement of the sites. The theory 
assumes that all sites have the same probability of forming bonds and that the correlation 
between adjacent sites is missing. 

The Helmholtz free energy of the system is written as a sum of the hard-sphere reference 
free energy, F HS , plus a bond contribution, F bond . The Helmholtz free energy due to bonding 
derives from a summation over certain classes of relevant graphs in the Mayer expansion^. In 
the sum, closed loops graphs are neglected. The fundamental assumption of the Wertheim 
theory is that the conditions of steric incompatibilities are satisfied: (i) no site can be 
engaged in more than one bond and (ii) no pair of particles can be double bonded. These 
steric incompatibilities are satisfied in both our models thanks to (i) the small 6 chosen for 
the short-ranged square-well attraction and to (ii) the location of the sticky sites on the 
hard-sphere particles surface. In the formulation of Ref.— , the bond free energy density of 
a system of /-functional particles is 

a rpbond i 

^-— = p ln(l-p b y + -pfp b (2) 

where (3 = p = N/V is the particle number density and p b is the bond probability. 

Since we assume equal reactivity for all sites, the bonding process can be seen as a chemical 
reaction between two unsaturated sites in equilibrium with a pair of bonded sites. In this 
respect one can write 

where T b is the site-site bond free-energy, i.e. the free energy difference between the bonded 
and the unbonded state. 

The Wertheim theory predicts an expression for T b in term of liquid state correlation 
functions and spherically averaged Mayer functions. More precisely 

where A refers to a single site-site interaction (since all bonding sites are identical) and it 
is defined as ^ 

A = 4vr f g H s(r 12 )(f(12)U^ 2 r 2 12 dr 12 . (5) 

Here gHs{i '12) is the reference hard-sphere fluid pair correlation function, the site-site Mayer 
function is /(12) = exp[— Vsw{ Y %) I^bT] — 1, and (/(12)) fa , 1)W2 represents an angle-average 
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over all orientations of particles 1 and 2 at fixed relative distance r\2- Since the Wertheim 
theory is insensitive to the location of the attractive sites, the number of interacting sites 
on each particle is encoded only in the factor / before A in Eq. |H For a site-site square-well 
interaction, the Mayer function can be calculated as 41 - 

(/(12)) wliWa = [exp(/M)) - 1] S(r) (6) 

where S(r) is the fraction of solid angle available to bonding when two particles are located 
at relative center-to-center distance r (r = 7*12), i.e. 



Sir) = g±£z£gg^£±I> . ( 7) 

ba z r 



The evaluation of A requires only an expression for gHs{ r ) i n the range where bonding 
occurs (a <r < a + 5). We use the linear approximation 42 

, s 1-0.50 9 0(1 + 0) /• - a 



9Hs{r) 



a 



(8) 



(1-0)3 2(1-0)3 

(where = |cr 3 p), which provides the correct Carnahan-Starling 4 ^ value at contact. This 
gives 



V h {eP U0 - 1} 

T 

( _ 5 (3a 2 + 85a + 3d 2 ) x 3 (\2ba + 55 2 ) 



A = ++"^x (9) 



2 a(15a + A5) Y 2 cr(15cr + A8) r 

where we have defined the spherically averaged bonding volume V& = 47r S(r)r 2 dr = 
7r(5 4 (15cr + 4<5)/30<t 2 . We note that the above expression of A simplifies in the low density 
limit. Indeed, when p — > 0, the hard-sphere pair correlation function tends to the ideal gas 
limit value gHs{ r ) ~ 1- in this limit A doesn't depend on p, i.e. A = Vb(e^ u ° — 1). We 
note that bonding takes approximatively place when exp(/3wo) ^ 1- Indeed bond formation 
arises from a balance between the energetic gain of forming a bond (A{7& = — mq) and an 
entropy loss (ASb), which is expressed in the theory as logarithm of the ratio between V& 
and the volume per bonding site, V / (fN)^. Since Vb <^V/ (fN), bonding becomes relevant 
when /3uq ^> 1. In the following we will thus approximate [e^ u ° — 1) with e^" to simplify 
the theoretical expressions. 

Once the free energy is known, it is possible to derive various equilibrium properties of the 
system through thermodynamic relations. We find expressions for the energy per particle, 



the specific heat maxima, the extent of polymerization and the pressure of the system in 
terms of pb, which is a function of T and p from Eq. [3j The potential energy per particle 
E/N is given by 

N = = -2 fUoPb (10) 

i.e. it is exactly the fraction of bond times —uof/2. The constant volume specific heat Cy 

can be calculated as 

n d (§ ) 1 f u 2 Q Pb(l-Pb) ( s 

Cv -^r~2 f f^ i+ Pb ■ 

At each p, the specific heat has a maximum (whose amplitude increases with /) at finite T, 
which defines a line of specific heat extrema in the (T, p) plane. The Cy ax line can be used 
as an estimate of the polymerization transition ii^Q^L-MM^dMS.. 

In the characterization of the self-assembly of particles, experimentalists often consider a 
quantity called extent of polymerization $(£), which is normally measured by spectroscopy. 
$ is defined as the fraction of particles bonded in clusters. This quantity plays the role of 
order parameter in the polymerization transition: it changes continuosly form the value zero 
at high T, when all particles are in the monomeric state, to the value one at low T, when 
particles are bonded in clusters. This crossover becomes sharper and sharper on decreasing 
p. Since the monomer density is simply obtained by the observation that all of the sticky 
spots on each particle must be unbonded, i.e. p\ = p(l — Pb)^ , the extent of polymerization 
is given by 

S = £^ = l-(l-p*)'. (12) 

As a function of pb, the branched polymerization transition becomes sharper and sharper 
on increasing the functionality of the system. 

The pressure P of the system can be evaluated by deriving, respect to the volume, the 
Wertheim free energy, i.e. (5P = —{df3F / 8V)t- The bonding contribution to P is thus 

"1 1 



P pbond 
=Pf 



2 l-p b _ 

In the low p limit (gHs( r ) ~ 1), it is possible to neglect the p dependence of A and (3P bond / p 
becomes equal to —\fpb- Appendix A provides a physical understanding of this expression. 
The hard-sphere contribution to the pressure can be evaluated via the Carnahan-Starling 
equation of stated 

PpHS (1 + + 2_ (/) 3 ) 

p (1-0)3 • 1 } 



From the resulting V and T dependence of P, it is possible to evaluate the liquid-gas 
coexistence region in the phase diagram, by solving the following set of equations 



T g = T l = T* 

P 9 = Pi = P* (15) 

fVg 

/ [P(V,T*) - P*]dV = 0, 
Jv l 

where T g , P g , V g and 7], Pi, V\ are respectively the temperature, the pressure and the volume 
of the two coexisting phases. The third equation corresponds to the Maxwell construction. 

The main assumption of the Wertheim theory is that molecules (or particles) cluster in 
open structures without closed bond loops. The hypothesis of absence of closed bonding 
loops is also at the heart of the Flory-Stockmayer theory, developed to model aggregation 
in chemical gelation. The Flory-Stockmayer theory 3 - provides expressions for the number 
density of clusters of n particles, p n , as a function of the bond probability (the extent of the 
reaction in the Flory-Stockmayer language). For functionality / 

Pn = p(i- Pb y [ptii-pty- 2 ] 71 - 1 ^ (i6) 

f(fn — n)\ 
Un ~ (fn-2n + 2)\n\ 

where p = ^2 n np n = N/V is the total number density. In Appendix A, we show that the 
Wertheim free-energy of Eq. (j2J) is equivalent to the free energy of a system of non-interacting 
clusters distributed according to Eq. [16j Here we make use of the Flory-Stockmayer theory 
for providing an expression, to be used in conjunction with the bond probability derived 
using the Wertheim theory, to evaluate the location in the (T, p) plane of the percolation 
line. The bond probability at percolation, is 

Pb = j^j- (17) 

Fig. [1] shows the resulting phase diagram evaluated according to the Wertheim theory 
for three different values of /. More specifically it shows the relative location of the phase 
coexistence line, the percolation and the maxima of specific heat line. According to the 
Wertheim theory, the coexistence region becomes wider on increasing /. For the case / = 5, 
at low T the gas coexists with a liquid with number density p w 0.8, a value significantly 
smaller than the one commonly observed for particles interacting via spherical potentials. 
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The percolation line merges into the coexistence curve on the left of the critical point, 
confirming that a spanning cluster of bond is a pre-requisite for the development of a critical 
phenomena^. For the shown / values, the locus of Cy ax is located below the corresponding 
percolation line, in agreement also with recent findings for a spherical model with / = 4M. 
However in the limit where / — > 2, realized via a mixture of / = 2 and / = 3 particles with 
average functionality 2.055—, the percolation line lies below the Cy ax line. The intersection 
of the Cy ax line with the coexistence curve progressively moves from the left to the right 
of the critical point on increasing /. Already for / = 5 the density at which the Cy ax line 
meets the coexistence line is more than twice the critical density. While it is not reasonable 
to extend the Wertheim theory to large / values, it is tempting to speculate that, on further 
increasing /, the intersection point will keep on moving to larger densities so that in the 
spherical limit case ( with analogous range of interaction) the entire Cy ax line lies in a 
physically inaccessible region (due to packing-driven kinetic arrest). 

IV. MONTE CARLO SIMULATION 

We perform simulations of the first model discussed in Sec. [Ill (in which the sticky spots 
location is regular) with the aim of evaluating the gas-liquid coexistence lines. We aim 
to provide a definitive proof that reducing valence generates a region of thermodynamic 
stability of the liquid phase down to vanishing temperatures in a wide range of densities. 
Previous studies of the same models were indeed focused only on the location of the critical 
points^. We perform Gibbs Ensemble Monte Carlo simulations (GEMC) in order to evaluate 
the phase coexistence region of one component systems with functionality /. The GEMC 
method^ 3 , allows us to study coexistence in the region where the gas-liquid free-energy barrier 
is sufficiently high to avoid crossing between the two phases. We simulate for about 5 million 
MC steps, where a MC step is defined as Aa = 10 5 attempts to translate and rotate a 
randomly chosen particle, Ajy = 10 3 attempts to swap a particle between the gas and the 
liquid boxes and Ay = 100 attempts to modify the volumes. A translational/rotational move 
is defined as a displacement in each direction of a random quantity distributed uniformly 
between ± 0.05 a and a rotation around a random axis of random angle distributed uniformly 
between ±0.1 radiant. The choice of such a large ratio between translation/rotation and 
swap attempts, Aa/Ajv = 100, is dictated by the necessity of ensuring a proper equilibration. 
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In the case of particles with short-range and highly directional interactions this choice is 
relevant, since the probability of inserting a particle with the correct orientation and position 
for bonding is significantly reduced as compared to the case of spherical interactions. 

We also study the model in which the sticky spots are non-regularly distributed on the 
surface (see Sec. |TT])- In particular we focus on the location of the critical point, since the 
values of critical temperature and density for the corresponding ordered arrangement have 
already been studied^. To assess the effect of the randomness on the location of the critical 
point we perform standard Grand Canonical (GCMC) simulations. In this ensemble, the 
chemical potential /i, the temperature T and the volume V are fixed. MC moves include 
insertion and deletion of particles as well as particle translation and rotations. Translational 
and rotational moves are identical to the one described above for GEMC. In each particle 
insertion move, a particle with a different realization of the location of the spots is placed 
in the box. GCMC simulations are extremely helpful in the study of the behavior of the 
system close to the critical point, since they allow for a correct exploration of the large 
range of densities experienced by systems in the vicinity of a critical point. To locate the 
critical point we perform simulations at fixed T, /i and V, and we tune T and /x until the 
simulated system shows ample density fluctuations, signaling the proximity to the critical 
point. Once a reasonable guess of the critical point in the (T, /z)-plane is reached, we start 
at least 8 independent GCMC simulations to improve the statistics of the fluctuations in the 
number of particles N in the box and of the potential energy E. The location of the critical 
point is performed through a fitting procedure associated to histogram reweighting 5 -^ and a 
comparison of the fluctuation distribution of the ordering operator Ai at the critical point 
with the universal distribution characterizing the Ising class 55 . The ordering operator Ai of 
the gas-liquid transition is a linear combination Ai ~ p + su, where p is the number density, 
u is the energy density of the system, and s is the field mixing parameter. Exactly at the 
critical point, fluctuations of Ai are found to follow the Ising model universal distribution 5 - 5 -. 



V. NUMERICAL RESULTS AND COMPARISON WITH WERTHEIM PREDIC- 
TIONS 

We first focus on the effect of patchiness on the phase coexistence region when the particles 
functionality / is small. Fig.[2]shows the numerical phase coexistence curves for systems with 
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a number / = 3,4,5 of attractive sites geometrically distributed on the particles surfaces. 
The figure clearly shows a strong reduction of the phase separation region, i.e. an extention 
of the region of stabilty of the liquid phase on decreasing /. Similar behavior is shown by 
the Wertheim predictions, despite the agreement gets worse on increasing /. Hence, both 
theory and simulations confirm 6,7 that a region of densities which is not affected by phase 
separation is a characteristic of patchy interacting particles systems. The reduction of the 
valence is thus crucial for suppressing the low temperature ubiquitous process of separation 
in a dense and dilute solution of particles always observed with spherical potentials. 

Fig. [2] also suggets that the small functionality of the particles makes it possible to observe 
chains and clusters in long-lived thermodynamic equilibrium. In other words, patchiness 
offers a way of sampling equilibrium homogeneous states in a large region of intermediate 
and small densities, where packing is not any longer the leading driving force controlling 
the structure of the system. The shrinking of the unstable region explains why particles 
interacting via a limited number of functional groups tend to form, at low temperature, 
open homogeneous structures, which are stabilized by an extended network of long-lived 
bonds. 

To assess if the shape of the coexistence does depend on / we show in Fig. [3] the same 
data of Fig. [2] represented as a function of reduced variables, T/T c and p/p c - While the 
Wertheim theory predicts a scaled width of the gas-liquid coexistence that shrinks with /, 
numerical data show that, far from the critical point, the curve for / = 3 appears to be 
significantly wider than the one for / = 4 and 5. Instead, close to the critical point the 
shape in reduced units appears to be rather insensitive on / (in agreement with previous 
findings 26 ) . 

Next we focus on the differences between a geometric and a random distribution of the 
patches and in particular on the / dependence of p c and T c in the two different cases. In the 
disordered case we vary / from 4 to 6. The results of the GCMC simulations are reported in 
Tab. [H together with corresponding quantities previously calculated for the geometric case 7 . 
The results for the critical point location in the two models are also graphically illustrated 
in Fig. |H The same trend with / is shown by both models. 

It is interesting to observe that, keeping / constant, T c and p c both decrease on moving 
from the geometric to the random arrangement of the sticky sites. This decrease suggests 
that the propagation of the connectivity is less efficient in the random patches case, speaking 
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for (i) a waste of bond formation possibilities and/or (ii) a failure in the development of long 
range paths of bonds. Concerning point (i), we note that a random distribution of patches on 
the particle surface may introduce correlation in the formation of adjacent bonds. Indeed the 
presence of a bonded interaction may induce a screening effect (and hence a decrease in the 
probability of forming bonds) on sites closeby located, due to excluded volume interactions. 
Concerning (ii), we note that a random distribution of sites may also favor the formation 
of closed loops of bonds, due to a increase in the number of angular possibilities which 
satisfy short ring structures, which are known to suppress the critical phenomenal. These 
observation can also explain why the Wertheim theory predictions (which are based on the 
assumption of both independent bonds and absence of ring structures) are closer to the 
geometric case model (see Tab. ITTI) . 

We note that Tab. [J and [Til also report the reduced values of the second virial coefficient 
at the critical point B^/B^ 3 . The analytical expression of B 2 / ' B^ s is the following 



where B 2 S = 2/3na 3 is the hard-sphere virial coefficient. 

We also evaluate the bond probability at the critical point, p%, on varing / in both 
the geometric and random patches models and we compare (see Tab. IIIII) the two sets of 
values with the Wertheim theoretical predictions. As previously observed, the Wertheim 
predictions are closer to the regular case, even if the theory is insensitive to the patches 
distribution on the particles surface. We also note that the critical bond probabilities in 
the geometric model are comparable with the ones recently calculated in Ref.— for particles 
with the same bonding geometry interacting via the Kern-Frenkel potential. On the other 
hand, yf h for the random model is significantly larger than for the ordered case, supporting 
our scenario of a less efficient propagation of connectivity in the random case as compared 
to the geometric one. 

Finally, we report in Fig. [5] the critical fluctuations distributions P(M) of the order 
parameter M. in both the geometric and random patches models with f — 5. The calculated 
distributions are compared to the expected fluctuations at the critical point for systems in 
the Ising universality class 5 -. The comparison provides evidence that the transition belongs 
to the Ising universality class in both studied cases. This is true for each studied value of /. 
The inset of Fig. [5] shows the corresponding density fluctuations distributions P(p) at the 



B 2 




(18) 
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estimated T c and critical chemical potential \i c . The distribution becomes more asymmetric 
on decreasing p c , signaling an increasing role of the mixing field s (see also Tab. [I]). This 
means that, at equal /, the density fluctuation distributions are more asymmetric in the 
random case rather than in the geometric one. 

VI. DISCUSSION AND CONCLUSIONS 

We study the /-dependence of the critical behavior in two different patchy models of 
/-functional particles. In both models, the patchy particles are hard-spheres decorated on 
their surface by a small number of identical sticky sites, interacting via a short-range square- 
well attraction. The difference between the two studied models is the arrangement of the 
attractive sites on the particles surface. In the first case sites are arranged on a regular 
structure (see Fig. 1 of Ref.~) in the same geometry of recently synthesized patchy colloidal 
particles^ 1 ^^. In the second case, the distribution of the sites is random and different for 
each particle. 

We compare numerical results and predictions of the thermodynamic perturbation the- 
ory developed by Wertheim2&>22. xhis theory assumes the condition of single-bond per patch 
and neglects the possibility of forming loops of bonded particles. As previously suggested 
in Ref.— , the free-energy expression provided by the Wertheim theory can be interpreted as 
the free energy of non-interacting clusters. We show in the Appendix A that the correspond- 
ing cluster size distribution is the one provided by Flory and Stockmayer in their seminal 
work on chemical gelation^^. In this respect, our study provides an effective expression 
for describing the density and temperature dependence of the free energy in self-assembly 
of branched structures and networks. The theory of equilibrium association for systems 
that form branched structures is receiving particular attention in the last yeaia^^^^^^, 
since these systems are found in many technological and biomedical applications, as well as 
in many biological processes. It is thus crucial to provide a general approach for describing 
the thermodynamics of the branched polymer self-assembly over the whole range of temper- 
atures, extending to branched system the work developed in the last decades for the case of 
self-assembling chains and wormlike micelles^ 1 ^ 1 ^. 

We explicitly solve the Wertheim theory for the chosen site-site interaction and we eval- 
uate lines of specific heat maxima (a signature of the presence of a specific bonding process) 
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and the the gas-liquid coexistence lines for 3 < / < 5. Thanks to the mapping between the 
Wertheim theory and the Flory-Stockmayer approach we also provide expressions for the 
dependence on / of the percolation line. We find that, for all studied /, the percolation line 
merges into the phase separation curve on the left hand side of the critical point, while the 
intersection between the Cy ax line and the coexistence curve moves from the left to the right 
of the critical point on increasing /. Even if the Wertheim basic assumptions can not be 
extended to high valence cases, we speculate that, on further increasing /, the intersection 
between coexistence and Cy ax line will further shift to larger densities. In this respect the 
absence of a Cy maximum in the spherical case could be due to the fact that the entire 
(jmax j- ne Y ies j n a re gion of large densities, made physically inaccessible by the progressive 
slowing down of the dynamics on approaching the glass transition. Indeed, on increasing /, 
the patchy potential tends to a spherical square-well model with analogous range of inter- 
action. For spherical potentials, it has been shown that the glass line — which provides the 
large-density limit of stability of the liquid state^ 1 ^ — intersects the coexistence line at a 
finite temperature and density. 

The Wertheim predictions for the the gas-liquid coexistence curve are compared to results 
of Gibbs Monte Carlo simulations of the regular patches arrangement model. We find that 
the reduction of the number of patches is accompanied by an enlargement of the region of 
stability of the liquid phase in the (T, p) plane, confirming the scenario suggested in Ref.-. 

Both the Wertheim theory and the simulations show that in models of reduced valence, 
states with uq ^> /c#T can be approached in equilibrium and reversibly. Thus in the pres- 
ence of patchy interactions it becomes possible in a wide range of densities to cool down 
the system progressively via a sequence of equilibrium homogeneous states. This is not pos- 
sible in spherically interacting particles for which phase separation always destabilizes the 
formation of a homogeneous arrested system at low T. Exploring homogeneous states at 
low temperature opens the way for sampling thermodynamic states characterized by bonds 
with very long lifetime. When the bond-lifetime becomes comparable to the experimental 
observation time, a dynamic arrest phenomenon at small packing fraction takes place. The 
reduction of the valence thus makes it feasible to approach dynamic arrest continuously 
from equilibrium and to generate a state of matter as close as possible to an ideal ge^. 
The relation between arrest in limited valence patchy colloidal particles and arrest in strong 
network forming molecular and atomic liquids have been recently discussed in Ref.— 
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Finally, we also study through Grand Canonical Monte Carlo simulations the location 
of the critical point for disordered arrangements of sites on the hard-sphere surface. Even 
in this case, we find that T c and p c moves towards lower temperatures and densities on 
decreasing the number of the patches. This confirms that the maximum number of bonds 
per particle plays an important role in controlling the stability of the liquid phase. The 
fact that the shift with valence of the critical point towards lower temperature and densities 
can be accomplished with either a geometric or random arrangement of patches could be 
particularly significant to experimentalists, since it indicates that ordered arrangements of 
patches are not absolutely necessary to achieve interesting assembly effects. 

We observe that, even if the Wertheim theory is insensitive to the arrangement of the 
sticky sites, the theoretical predictions for the critical point location in the phase diagram 
are closer to the geometric case model rather than the random one, suggesting that the 
propagation of the connectivity is less efficient in the random patches case. We recall that the 
theory is based on the assumption of (i) independent bonds and (ii) absence of closed loops 
of particles. Hence, the reduced connectivity of the random model, at equal temperature 
and density, could be related to (i) a reduction of bond formation possibilities, induced by a 
correlation between nearby sites, and/or to (ii) an increase in possibilities of ring strucures 
formation, which disfavor the development of branched bonding patterns. 

As a side remark, we add in Appendix B the demonstration that, within the Wertheim 
theory, when particles interact only via bonds and no hard-core repulsion is present, the 
thermodynamic stability line (spinodal) coincides with the percolation line. 
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VII. APPENDIX A 

Here we provide a physical insight of the Wertheim theory, by discussing two equivalent 
alternative derivations of the Wertheim bond free-energy. Both derivations assume the 
system of associating particles to be formally equivalent to a system of non-interacting 
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clusters, in thermodynamic equilibrium. For simplicity, we assume a De Broglie length 
A = a — 1 in both derivations. 

In the first derivation, we assume that the cluster-size distribution is provided by the 
Flory-Stockmayer expressions (see Eq. fTBl) . i.e. that the system of N /-functional associating 
particles aggregates in clusters characterized by the absence of closed bonding loops. Bonds 
are also assumed to be uncorrelated so that to each bond is associated the same single-bond 
free energy Tb- In the absence of loops, the number of bonds in a cluster of size n is exactly 
(n — 1), since each new bond adds one new particle. Hence the bond free energy of the 
cluster is (n — If clusters do not interact, the system free energy F can be written as 

sum of the ideal-gas free energy of each distinct bonding topology cluster-type (accounting 
for the translational center of mass degrees of freedom) and a sum over the cluster bond 
free energies. Defining p^ as the number density of clusters with size n and with bonding 
topology k one can write 

y = EE^M-i]+EE p k n( n - wn- ( 19 ) 

n k n k 

Here V is the volume, (3 = 1/ksT (with k^ the Boltzmann constant) and the sum on n runs 
over all the possible cluster sizes, from one (monomers) to infinity, while the sum over k 
includes all u) n distinct bonding topology of clusters of size n. Since clusters with the same 
size but different bonding pattern are equiprobable, then p\ = p n /u) n (see Eq. [16]). Thus 
Eq. [19] becomes 

^ = VV [h^-ll +y> n (7i-l)#F6. (20) 

Substituting Eq. [16] in Eq. [201 and summing over n one obtains, for £>&<(/ — 1) _1 (which 
express the condition that all clusters are finite 39 ), the following expression for the system 
free energy in term of pb and bond free energy 

3F 

y = plnp-p+ (21) 



+ pln(l -p b 



— n, \ J 



ln P(l-^) e -^_ 1 
Pb 



2 

Such expression can be seen as a high temperature contribution^ (plnp — p) plus a remain- 
ing bonding term. The bonding free energy coincides with the Wertheim expression^^lL^ 
(pln(l —pbY + ^rp-) when the connection between pb and Tb is given by Eq. [3] This simple 
derivation can be also extended to binary mixtures. 
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An even simpler derivation has been suggested in Ref.— and it is here reported for com- 
pleteness. Also this derivation assumes that the system is an ideal gas of clusters and hence 
that the product (3PV is identical to the number of clusters N c . The evaluation of N c is 
straightforward for p b values smaller than the percolation threshold pf, since, in the absence 
of closed bond loops, N c is equal to the number of particles minus the number of bonds N b . 
Calling N™ ax = ^ the maximum number of possible bonds and noting that p b is the ratio 
between N b and N™^ , one finds 

N c = N - N b = N (l - Lp h ) (22) 



pp = p[l-lp b ). (23) 



and 

BP = n (l - ^ 

Since the system is in dynamic equilibrium, the particle chemical potential is independent 
from the cluster to which the particle belongs to. Hence, it is identical to the chemical 
potential \x of the monomer. The ideal-gas hypothesis implies that the activity of the 
monomer z = exp(/?/i) is related to the monomer number density by z = p\ = p(l — p&) ■ 
Hence the system free energy density can be immediately written as 

(3F = p(3p-(3P = p In [p(l - Pb )f] -p(l- f -pA (24) 

which coincides with the Wertheim expression when the reference system is the ideal gas. 

We note in conclusions that the above relations are valid only before percolation. Indeed, 
the sums over n in Eq. [20] as well as Eq. [22] are valid only for p b < p p b . Hence, the region 
of validity of the Wertheim theory should be strictly limited to non-percolating states. 
Nevertheless we observe that the theory works well even below the percolation threshold^ 1 ^. 
It could be in principle possible to extend the formalism to p b > pi by accounting correctly 
for the p b dependence of the number of clusters (which is always possible, analytically for 
small / and numerically above), but it is not clear how to handle the free energy contribution 
of the percolating cluster. 

VIII. APPENDIX B 

In this Appendix we examine the thermodynamic stability of a system composed of non- 
interacting clusters, described by the free energy of Eq. [20] A stable system is characterized 
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by a negative volume derivative of the pressure. The region of stability is delimited by the 
so-called spinodal line, defined as the locus of points such that (d/3P/dV) T = 0. The volume 
derivative of the pressure, under the ideal gas approximation gHs{ r ) = 1, is controlled only 
by the volume derivative of pi, (see Eq. 125]) . Interestingly, it gives for the bond probability 
at the spinodal line pi 

i.e. the same condition that defines percolation. Hence p s b = p\. The system is thus 
mechanically stable only in the non-percolating region. In other words in the ideal gas 
approximation no dense stable states are possible and the system exists only in the gas- 
phase. In the Wertheim theory the existence of a liquid phase is generated by the significant 
increase of the pressure at low V introduced by the hard-sphere reference contribution. Fig. O 
provides an example of the effect of the hard-sphere contribution on the pressure for the 
case / = 3. 

We conclude noting that the absence of the hard-core repulsion appears to be essential in 
formally associating the percolation line with the spinodal line, providing an analytic simple 
example of the possibility of interpreting critical phenomena in term of percolation^^. 
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/ 


T c 


Pc 


Pc 


s 


L 




3 geometric 


0.094 


0.141 


-0.471 


0.46 


9 


-28.772 


4 geometric 


0.118 


0.273 


-0.418 


0.08 


7 


-5.080 


5 geometric 


0.132 


0.351 


-0.410 





7 


-2.866 


4 random 


0.102 


0.208 


-0.531 


0.46 


8 


-21.978 


5 random 


0.118 


0.258 


-0.500 


0.25 


8 


-8.500 


6 random 


0.133 


0.310 


-0.482 


0.22 


7 


-4.258 



TABLE I: Values of the relevant parameters at the critical point for the geometric 7 (/ = 3, 4, 5) 
and random (/ = 4, 5, 6) cases: T c is the critical temperature, p c is the density of the critical point, 
(jl c is the critical chemical potential, s is the field mixing parameter and L indicates the largest 
studied box size. B^/BZ? 3 is the value of the reduced second virial coefficient at the critical point. 
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/ 


T 


Pc 




3 theory 


0.0925 


0.086 


-34.378 


4 theory 


0.1121 


0.154 


-8.498 


5 theory 


0.1275 


0.212 


-4.052 


6 theory 


0.1411 


0.261 


-2.414 



TABLE II: Critical values of the temperature and density for 3 < / < 6 evaluated through the 
Wertheim theory. We also report the corresponding values of the reduced second virial coefficient. 
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/ 


pi theory 


p^ geometric 


p\ random 


3 


0.633 


0.728 




4 


0.488 


0.640 


0.737 


5 


0.417 


0.577 


0.615 


6 


0.360 




0.539 



TABLE III: Critical values of the bond probability pb for / varying from 3 to 6. Theoretical values 
are obtained solving Eq. [3] at the critical point, while numerical ones are obtained as the ratio 
between the potential energy at the critical point and the energy of the fully bonded system. 
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FIG. 5: 
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FIG. 6: 
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FIG. 1: Theoretical predictions for the phase diagram of patchy systems on varying the particles 
functionality form / = 3 to 5. Coexistence curves and Cy ax lines are evaluated according to the 
Wertheim theory, respectively from Eq. [JS] and by finding the zeroes of the temperature derivative 
of Eq, UH i.e. [dCy / dT)y = 0. Percolation lines are evaluated according to the Flory-Stockmayer 
theory as the locus of points in the (T, p) plane such that pb(T, p) = j? h (see Eq. [TT|h 

FIG. 2: Gas- liquid coexistence regions in the (T, p) plane on varing / from 3 to 5. Points are 
GEMC numerical results for the model in which the sticky sites are geometrically arranged on 
the particles surface. Solid lines are Wertheim theoretical predictions for the coexistence curves 
obtained from Eq. [T5j The numerical (stars) and theoretical (crosses) critical points frorn^ are 
reported to help visualizing the critical behavior. 

FIG. 3: Gas-liquid coexistence curves in terms of the reduced temperature T/T c and the reduced 
density p/p c for systems with / = 3,4, 5 sticky sites. Points are GEMC results for the geometric 
arrangements of the patches while solid lines are the Wertheim theoretical predictions. 

FIG. 4: Comparison between numerical results for patchy particles with different number of sticky 
spots per particle in a geometric (squares) and random (circles) arrangement. Panel (a) shows the 
location of the points in the (T, p) plane. Panels (b) and (c) compare respectively the / dependence 
for T c and p c . Data for the ordered case are reproduced from Ref. 7 . 
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FIG. 5: Comparison between the critical fluctuations distributions of Ai ~ p + su in both the 
geometric and random cases with / = 5. The calculated P{M) are compared to the expected 
distributions (full line) for systems at the critical point belonging to the Ising universality class^. 
The inset shows the comparison between the corresponding density fluctuations distributions P (p) 
in the two cases. 

FIG. 6: Equation of state for two different isothermes (T = 0.08 and 0.1) in the pressure vs. 
volume plane in the case of three functional particles. Solid lines (id) are related to the ideal 
system described by Eq. [23] Note that in the high T limit Eq. [23] reduces to the ideal gas equation 
of state {pPjp = 1). Points on solid lines indicate the maxima of the pressure respect to p -1 : the 
set of these points on varing T provides the spinodal line of the system. Dotted lines (hs) are related 
to a system in which bonding is accompained by an hard-core repulsion: the bonding contribution 
is given by Eq. [13] in which the density dependence of the <?j/s(r) (see Eq. [8]) is considered, while 
the hard-sphere contribution is given by Eq. [T^l 
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